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I. INTRODUCTION 

Usually the dominant interatomic interactions in an 
atomic Bose-Einstein condensate (BEC) are asymptoti- 
cally of the van der Waals (vdW) type, which falls off as 
r -6 and is short-range in comparison to the de Broglie 
wavelength of the atoms. These interactions can be incor- 
porated into a mean-field description of the condensate 
via a delta- function pseudo-potential P, 

(1) 



U(r) = 4irh 2 a s S(r)/m = gS(r) 



involving just the s-wave scattering length a s and atomic 
mass m. The interactions then appear as a cubic non- 
linearity in the Gross-Pitaevskii equation [3j for the order 
parameter yb(r) of a trapped zero-temperature BEC 



Mr) = \ 2^ V +Wr)+<#(r)| W) > ( 2 ) 

where /i is the chemical potential. The trapping poten- 
tial, Vt ra p, due to a magnetic or optical trap is typically 
harmonic, Vt ra p = (m/2)[w 2 x 2 + to 2 y 2 + uj 2 z 2 }. 

The Thomas-Fermi regime for a trapped BEC is 
reached when the zero-point kinetic energy is vanishingly 
small in comparison to both the potential energy due to 
the trap and the interaction energy between atoms 0- 
Many of the current atomic BEC experiments satisfy 
these conditions, which tend to hold for condensates con- 
taining a large number of atoms. When the kinetic en- 
ergy is neglected the time-independent Gross-Pitaevskii 
equation J5J can be trivially solved for the static conden- 
sate density, 

n(r) = |V(r)| 2 = [p - V tlap (r)]/g for n(r) > 0, (3) 

and n(r) = elsewhere. Thus the density profile is com- 
pletely determined by the trapping potential, and in a 
harmonic trap n(r) has an (inverted) parabolic profile 
and the same aspect ratio as the trap. 

Our aim here is to obtain similar exact results for a 
BEC in which dipole-dipole interactions play an impor- 
tant role. Compared to the vdW interaction the dipo- 
lar interaction is long-range and anisotropic, and conse- 
quently these systems can be expected to exhibit novel 



behavior including unusual stability properties _ 
exotic ground states such as supersolid 0,13 and checker- 
board phases , and modified excitation spectra 0, El , 
even to the extent of a roton minimum |llL Il2| in the dis- 
persion relation. We have recently reported on the exact 
dynamics of a dipolar BEC in the Thomas-Fermi limit 
[l3j , including the quadrupole and monopole shape oscil- 
lation frequencies. Here we give the full derivation of the 
static solution, which was only stated in |13| , and investi- 
gate stability and electrostriction (change in volume due 
to an applied electric field) . We also calculate the dipole- 
dipole potential outside the boundary of the condensate, 
which has a bearing on the distribution of thermal atoms 
and on the stability of a dipolar BEC, but also on a lat- 
tice array of dipolar BECs, since the effective giant dipole 
on each site is coupled to its neighbors by dipole-dipole 
interactions [T3 |. 

The long-range part of the interaction between two 
dipoles separated by r, and aligned by an external field 
along a unit vector e, is given by 



E/ d d(r) 



C 



dd , „ (Sjj 

4ir eiBj 



(4) 



where the coupling Cdd/(47r) depends on the specific re- 
alization. Marinescu and You investigated the low- 
energy scattering of two atoms with dipole-dipole inter- 
actions induced by a static electric field, E = Ee, so that 
Cdd = E 2 a 2 /eQ. Yi and You then went on to con- 
sider a BEC composed of such atoms. Another possible 
scenario is permanent magnetic dipoles, d m , aligned by 
an external magnetic field, B = Be, leading to a coupling 
Cdd = Mo^m- A BEC with magnetic dipole-dipole inter- 
actions was first discussed by Goral, Rzazewski, and Pfau 
J4J. A measure of the strength of the long-range dipole- 
dipole interaction relative to the s-wave scattering energy 
is provided by the dimensionless quantity 



c 



dd 



3.9 



(5) 



This definition arises naturally from an analysis of the 
frequencies of collective excitations (Bogoliubov spec- 
trum) in a homogeneous dipolar BEC |4.ll7|. As we shall 



2 



see later on, the BEC is stable as long as < £dd < 1, 
but loses that stability in the Thomas-Fermi limit when 
£dd > 1- Note that in the presence of a strong electric 
field the s-wave scattering length can be modified |sL fl~5l| . 
and therefore g and hence £dd should be treated as ef- 
fective quantities when dealing with electrically induced 
dipoles. 

Although the magnetic interaction between two atoms 
is often masked by a stronger s-wave interaction, two re- 
cent proposals indicate how magnetic dipolar effects can 
be enhanced, either by rotating the magnetic field in res- 
onance with a collective excitation frequency of the sys- 
tem 01 > or by using a Feshbach resonance to reduce the 
s-wave scattering length [Isj . Other suggestions to real- 
ize BECs with strong dipolar interactions include polar 
molecules |4| and Rydberg atoms |5| • Laser induced (dy- 
namic) dipole-dipole interactions differ from the static 
case of Eq. Q by extra retarded terms, including an 
term which is always attractive; they are discussed in 
References and [2£j. In certain situations this very 
long-range part of the interaction is important and can 
be responsible for unique features such as self-binding, 
and plasmon-like collective excitations. Here, however, 
we confine ourselves to the static case, which could be 
realized by static fields but also by using a laser provided 
the atomic separation is considerably smaller than the 
laser wavelength. 



II. THOMAS-FERMI EQUATION FOR A 
DIPOLAR BEC 

Proceeding from the Thomas-Fermi equation for a 
static BEC, i.e. the time-independent Gross-Pitaevskii 
equation without the kinetic energy term, we seek an ex- 
act solution for the density, n(r), of a condensate with 
both dipole-dipole interactions and the usual short-range 
s-wave scattering in a harmonic trap, which for simplicity 
(though not necessity) we take to be cylindrically sym- 
metric {u) x — uj y ). The Thomas- Fermi equation then 
reads 

P=\ m Hp 2 + ^) + 9<r) + $ dd (r) (6) 

where p 2 — x 2 +y 2 , and $dd(r) is the mean- field potential 
due to dipole-dipole interactions 



$ dd (r) = J dV U dd (r - r')n(r') 



(7) 



The intuitive form of Eq. Q, which corresponds to the 
Born approximation for two-body scattering, has been 
shown by Yi and You to be accurate for dipole mo- 
ments of the order of a Bohr magneton and collisions 
away from any shape resonances. Recently Derevianko 
PH proposed a more sophisticated approach to the dipo- 
lar scattering problem which suggests that dipole-dipole 
interactions can be substantially larger than previously 
estimated 22] . However, it appears that, provided the 



condensate is not too strongly deformed, the basic form 
of Eq. |(7J| with the bare interaction Q remains valid, 
albeit with a renormalized coupling Cdd H3 

The presence of the non-local dipolar mean-field poten- 
tial <i>dd(r) means that the Thomas- Fermi Eq. JfjJ) is an 
integral equation and so less trivial to solve than in the 
purely local delta-function pseudo-potential case. How- 
ever, it is straightforward to demonstrate that this equa- 
tion also admits an inverted-parabola as a self-consistent 
solution. We begin our analysis with a suggestive re- 
casting of the dipole-dipole term using the mathematical 
identity 



[6ij - 3fifj) 



We can then write 



1 4tt 

r 



-ViV ir --3-MW- (8) 



<KMv) = Cm <v, ( V<V^(r) + ^n(r) ) 



with 0(r) = — f 
4vr J 



1 f d r' n(r' 



(10) 



The problem thereby reduces to an analogy with electro- 
statics, and we need only calculate the 'potential' </>(r) 
arising from the 'static charge' distribution n(r). In par- 
ticular, <j)(r) given by (|10|l must obey Poisson's equation 
V 2 = — n(r). We adopt the following inverted-parabola 
as an ansatz for the density profile of an iV-atom conden- 
sate 



n(r) = n 



i _ £. _ il 

r 2 m 



for n(r) > (11) 



with radii R x = R y and i? z , and where the central density 
no is constrained by normalization to be 



n = 15N/(8ttRIR z ) . 



(12) 



Then Poisson's equation is satisfied by an 'electrostatic 
potential' of the form 

0(r) = a a + a±p 2 + a 2 z 2 + a a p 4 + a^z A + a 5 p 2 z 2 . (13) 

However, by Eq. (jSJ), the physical dipolar contribution 
^ddlr) to the mean-field potential inside the inverted- 
parabola BEC (fill) will now itself also be parabolic, just 
like the potentials due to the harmonic trap and the local 
s-wave scattering interaction. Thus the Thomas-Fermi 
equation |J?jJ contains only parabolic and constant terms 
and so, remarkably, just as in the pure s-wave case, in 
the presence of dipole-dipole interactions a parabolic den- 
sity profile is also an exact solution of the Thomas-Fermi 
problem in a harmonic trap, although this time we should 
expect that the condensate aspect ratio differs from that 
of the trap. It remains to determine the coefficients ap- 
pearing in (| 1 31) and in (|11|) and adjust them in such a way 
that the Thomas-Fermi equation is satisfied. To this end 
we shall evaluate the integral (|10|) for a density of the 
form (|ll|l . This is an arduous task because the domain 



of integration is bounded by and has the symmetry of a 
spheroid or, in the general case, even of an ellipsoid. Cal- 
culating the integral is possible only if one takes explicit 
account of this symmetry, and we shall demonstrate two 
independent ways of doing that. One is to transform into 
spheroidal coordinates, use the known Green's function 
of Poisson's equation in these coordinates [23 , and subse- 
quently transform back into Cartesian coordinates. The 
other is to start from basics and integrate over successive 
thin ellipsoidal shells. While the former approach is also 
quite involved it is simpler than the latter. However, if 
we were to drop our simplifying assumption of cylindrical 
symmetry, the second approach is the only workable as 
the general solution of Poisson's equation in general ellip- 
soidal coordinates is unmanageably complicated because 
the se par ation constants do not separate in these coordi- 
nates |23l p. 757]. The second approach is presented in 
appendix \K\ 



III. GREEN'S FUNCTION IN SPHEROIDAL 
COORDINATES 

We now demonstrate the Green's function approach. 
For prolate spheroidal coordinates (£, 77, ip) we have x = 

q\/(£ 2 ~ !)(! - V 2 )costp, y = qy/(£ 2 - 1)(1 - r? 2 )sin^, 
z = q£r/. Surfaces of constant £ are confocal spheroids 
whose eccentricity is l/£, and £ runs between 1 and 
00. Surfaces of constant 77 are confocal two-sheet hy- 
perboloids of revolution, and 77 runs between -1 and 1. 
For R z > R x the boundary of the density profile 
is a prolate (cigar-like) spheroid with semimajor axis 
R z , semiminor axis R x , and eccentricity yl — R 2 /R 2 . 
To make the spheroidal coordinate system confocal to 
that boundary we need to choose the scaling constant 
q = y/ ' R 2 — R 2 . Then we can use the Green's function in 
prolate spheroidal coordinates [2j| to write the potential 
EH as 



p2 p2 



d£' /* V(£' a -» 7 ,2 )n(£W) 
1 J -1 



^=0.8 



^=0.4 




T|=COS Jl/3 









X 




T|=-COS 7l/3 



FIG. 1: An illustration of oblate spheroidal coordinates. Sur- 
faces of constant £ are confocal spheroids, with £ = being 
a flat disk of radius q and £ — > 00 an infinitely large sphere. 
Surfaces of constant \rj\ are one-sheet hyperboloids (narrow- 
waisted reels); 77 changes sign with z. The surface described 
by 77 = is the xy plane with a circular hole of radius q, and 
\n\ = 1 is the z axis. 



r 2 )/{2q) and 77 = (n. - r 2 )/(2q) with n = [x 2 + y 2 + 
(z + q) 2 ] 1 / 2 and r 2 = [x 2 + y 2 + (z — q) 2 ] 1 ^ 2 - We thereby 
obtain a 'potential' of the form predicted by Eq. I|13l) : 



(r)= "°^ 9 J 245(l-„ 2 ) 2 
V ' 192(1 — /v ) V 



+48(l-^)(2-H)( — 



-24(1-k 2 )(2-k 2 S) 



2^ I _P_ 
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3(2« 2 ~8 + 32)( — 



-3[2(2-5 K 2 ) + 3 K 4 S] [-£- 



+24(2 + An 2 - 3k 2 S) 



, 2 / \ 2 
p \ I z 



R. r 



R* 



(14) 



de / d V '(e 2 -v ,2 )n(e,v') 



where K = R x / R z is the aspect ratio of the BEC and 



1 



vT 



In 



1 + VT^ 



1-VT 



where Pi are Legendre functions of the first and Qi of 
the second kind. Since n(r) is quadratic in x, y, and z 
it is quadratic in £ and 77, and all integrals in the above 
expression are elementary. Performing the 77' integration 
first we see that the only contributing I are 0, 2, and 
4. To re-express the result for </>(£, 77, (p) in Cartesian 
coordinates we need to make the substitutions £ = (ri + 



for k < 1 (prolate). 

(15) 

If R x > R z then the boundary of the density 
profile l|ll|) is an oblate (pancake-like) spheroid, and 
we have to use oblate spheroidal coordinates x = 
Q\/(£ 2 + 1 )( 1 - 7 ? 2 ) cos V, V = 9v / (C 2 + 1 )( 1 -?7 2 )sm^, 
z = qS, 1 !- Surfaces of constant £ are again confocal 
spheroids but now with eccentricity 1 / + 1, and £ 
running between and 00. An illustration of oblate 
spheroidal coo rdinates is given in Fig. ^ We have to 
choose q = \J R 2 — R 2 to make the coordinate system 



confocal to the boundary of n(r). Using the Green's func- 
tion in oblate spheroidal coordinates |23l l24j we find for 
the potential 



xi„ — it, 



V t dr/(£' 2 +7/ 2 M£W) 



d£'/ cW 2 + 7/ 2 )ra(£W 



i^(2£+l)^( 77 )P,(7 7 ')^(iC)QH^') 



To return to Cartesian coordinates we 
need to make the substitutions £ = 
(v/x 2 + y 2 + (z + iq) 2 + x 2 + y 2 + (z - iq) 2 )/(2q) and 
r] = (^x 2 +y 2 + (z + iq) 2 ~^x 2 + y 2 + (z - iq) 2 )/(2iq). 
Then we find that the result for the potential is the 
same as in Eq. H14JI but with 
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FIG. 2: Aspect ratio of the condensate as a function of the 
dipole-dipole to s-wave coupling ratio £dd- Each line is for a 
different trap aspect ratio 7, which can be read off by noting 
that n{edd = 0) = 7. When < k < 1 the condensate is 
prolate, for k > 1 it is oblate. Likewise, when < 7 < 1 the 
trap is prolate, and when 7 > 1 the trap is oblate. Dashed 
lines indicate unstable branches. 
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arctan \/ k 2 — 1 for k > 1 (oblate) 



(16) 

The prolate and oblate cases are of course connected by 
analytic continuation, which however cannot be used to 
determine one from the other because of the ambiguity 
of sheets in the complex plane. 

In order to simplify the expressions that will follow, 
and in order to conform with existing notation in the 
literature 0, l25| , rather than working with the function 
S(k), we shall work instead with f(n): 



/(«) 



2 + k 2 [4-3S(k)] 
2(1 - k 2 ) ' 



(17) 



/(k) is a monotonically decreasing function of k, having 
values in the range 1 > /(k) > —2, passing through zero 
at k = 1. 



IV. SOLUTION OF THE THOMAS-FERMI 
EQUATION 

In this paper we shall take the external field to be along 
the z-axis. Then the result (|14fl for the 'electrostatic 
potential' </>(r) yields, by virtue of Eq. ©, a parabolic 
dipolar potential 



$dd = 



n C, 



dd 



Rl R 2 J(K> \ 2R1- R 2 



(18) 



which is valid inside the condensate (the potential outside 
the condensate boundary will be discussed in the next 
section). Substituting $dd(r) into the Thomas- Fermi Eq. 
© and comparing the coefficients of p 2 , z 2 , and 1, yields 



three coupled equations. The first equation, due to the 
constant terms, gives the chemical potential 



/1 = gn [1 - £dd/(«0] 



(19) 



This equation indicates that, all other things being equal, 
the effect of dipole-dipole interactions is to lower the 
chemical potential (which is proportional to the mean- 
field energy per particle) of a prolate (k < 1) condensate, 
whilst raising that of an oblate (k > 1) condensate. The 
radii R x (= R y ) and R z of the exact parabolic solution 
(|llfl are obtained from the coefficients of p 2 and z 2 . Wc 
find 



Rx = Ry = 



IbgNn 
Aitmuj 2 



1 + £dd 



2 1 — K 2 



— 1 



1/5 



and R z = R x /n. The aspect ratio n is 
solving a transcendental equation 

3K 2 £dd ' ' 1 



(20) 
determined by 



1-K 2 



+ ( £dd - 1) (n 2 - 7 2 ) = 
(21) 

where 7 = uj z /uj x is the ratio of the harmonic trapping 
frequencies. In fact, a property such as the aspect ra- 
tio is insensitive to the details of the density profile and 
Eq. (|21|l has been obtained previously from a Gaussian 
variational ansatz for the density |25| ■ Figures |21 and 
13 show examples of the dependence of k upon £dd for 
oblate, spherical and prolate traps. The effect of dipole- 
dipole forces polarized along the z-axis is to make the 
condensate more cigar-shaped along z. For an oblate 
trap (7 > 1) the BEC becomes exactly spherical when 
e dd = (5/2)( 7 2 -l)/(7 2 + 2). 

In order to illustrate the static properties of the 
Thomas-Fermi solution for a dipolar BEC we imagine an 
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FIG. 3: Aspect ratio of the condensate as a function of eaa- 
This is an expanded version of Fig. [5] illustrating the sur- 
prising result that beyond a certain critical oblateness of the 
trap (7 > 7 C rit = 5.1701) the system is metastable to scaling 
deformations, even for arbitrarily high values of the dipolar 
interaction strength (but is not necessarily stable to other 
types of perturbations such as phonons — see later). At the 
boundary value 7 = 7 C rtt, one finds k — > 2.5501 for £dd — > 00. 



experiment with a large fixed number of atoms, N, in a 
trap set to a particular aspect ratio, 7, where the value of 
£dd is adiabatically increased from zero. For electrically 
induced dipoles this would involve increasing the electric 
field, whereas for magnetic dipoles one could either rotate 
the external magnetic field, gradually changing the angle 
of rotation |17j . or reduce the s-wave scattering length 
using a Feshbach resonance. The system then follows 
one of the curves shown in Figures |2] and |21 In the ab- 
sence of the external field, when Edd = 0, the condensate 
aspect ratio matches that of the trap, k = 7. When the 
dipole-dipole interactions are switched on the condensate 
becomes more prolate than the trap and one always has 
K < 7. As long as < e^d < 1, the transcendental equa- 
tion (|21[l has a single solution, k, for any choice of trap, 
7. The behavior for Edd > 1 is more complicated and re- 
quires an analysis of the stability properties of a dipolar 
BEC, which we give in the next but one section. 



V. DIPOLAR POTENTIAL OUTSIDE THE BEC 

For a variety of applications it is very useful to know 
the potential outside a dipolar condensate. For example, 
in an array of dipolar BECs on a lattice the condensate 
at each lattice site can behave as a single mesoscopic 
spin [lj|. In order to determine the spin-spin coupling 
between sites one needs to know the external potential 
generated by each condensate. In Section lVl Cl below we 
shall see that the knowledge of the outside potential also 
gives insight into the problem of the stability of a single 
dipolar condensate. 

In order to calculate $dd(r) in Eq. (7J outside the con- 
densate, we again use relation JSJl and write the outside 



dipole-dipole potential in the same way as in Eq. ©, ex- 
cept that the term with does not arise because n(r) is 
obviously zero outside the condensate. Using the Green's 
function in prolate spheroidal coordinates we find for the 
potential (|10|) outside 



E>2 p2 

xt_ — -ft... 



d£7 W 2 W 2 ) 



OO 



where, as before, only i = 0,2,4 actually contribute to 
the sum. In the oblate case a similar formula applies, 
with just £ and £' replaced by i£ and i£' and the £' inte- 
gration running from to 1/Vk 2 - 1. We obtain 



(outside) 
dd 



(r) 



3cf£dd«o>« 2 



4 ( 1 _ K 2)3/2 

+ [9£V- 3(£ 2 + J7 2 ) + 
in the prolate case, and 

(outside) 3g£ddHoK 2 



6£(1 - 3/7 2 



(22) 



l] In 



<i<i 



(r) = - 



{6C(1- 3T7 2 ) 



4 ( K 2 _ 2)3/2 

[9£V - 3(£ 2 - if) - 1] (tt - 2 arctan^)} 



(23) 
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ily 



the oblate case. These expressions can eas- 
be converted back from prolate or oblate 
spheroidal coordinates (£, 77) into Cartesian coordi- 
nates, as was described above, by substituting £ 



V 



+ (z + q) 2 + 



.'/- 



[z_ 

c 2 + y 2 
coordinates, 



y* + (z 



g) 2 )/(2g ) and 

9) 2 )/(2<z) 
and £ = 



T) = {\Jx 2 + y 2 + (z + q) 2 
for prolate spheroidal 
{\/x 2 +y 2 + (z + iq) 2 + v / x 2 +y 2 + (z - iq) 2 )/{2q) and 
V = Wx 2 +y 2 + {z + iq) 2 - ^x 2 +y 2 + {z-iq) 2 )/{2iq) 
for oblate spheroidal coordinates. 

While easy to obtain, the expression for $^° utslde ) ( r ) 
is rather lengthy and unwieldy in Cartesian coordi- 
nates. Thus, if one has to work in Cartesian coor- 
dinates, one may prefer the asymptotic expression for 
the dipole-dipole interaction potential at large distances 
r = (p 2 + z 2 ) 1 / 2 , which is approximately 

N \( 1 _3z^ 
Anr 3 \ r 2 



(outside) 



dd 



(r) ~ Cdd 



R2 7j2 



45 z 



15 z A 



i4-y^ + y- 1+0 



R x , R z 



(24) 



and holds for both prolate and oblate condensates. It 
turns out that this asymptotic expression serves remark- 
ably well even quite close to the condensate. We give 
a derivation from integration over thin ellipsoidal shells 
in Appendix B. Note that equation l|24|l says that to a 
first approximation, when seen from outside the dipolar 
condensate behaves like a single giant dipole of N-times 
the single-atom dipole magnitude. The higher multipoles 
depend on the shape of the BEC. 
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VI. STABILITY OF A DIPOLAR BEC 

The partially attractive nature of the dipole-dipole in- 
teraction has been widely predicted 0, IE IE 13 to 
lead to a collapse of the BEC when the dipolar interac- 
tion strength exceeds a certain critical value, e™^. In the 
Thomas-Fermi limit depends only upon the trap as- 
pect ratio 7. Both the parabolic solution presented here 
and the Gaussian variational ansatz indicate that above 
the system is liable to collapse towards an infinitely 
thin and long prolate 'pencil' oriented along the field po- 
larization direction i.e. n — > 0, since the system lowers its 
energy by arranging the dipoles end-to-end. However, in 
reality a transition to another (more structured?) state 
(see H S IH 113) presumably occurs in preference 
to the system becoming truly singular. Bearing this in 
mind we shall consider below three nominally different 
types of instability: local density perturbations, 'scaling' 
deformations, and the 'Saturn-ring' instability. The lat- 
ter occurs due to a peculiarity in the potential that would 
be seen by an atom located outside the boundary of the 
BEC, and may result in a previously unforeseen type of 
instability. All of them predict the onset of instability 
when £dd > 1- Nevertheless, we have included in the fig- 
ures values of £dd exceeding unity, our justification being 
partly mathematical curiosity, and partly the fact that 
the inclusion of kinetic energy would extend the stability 
of a dipolar BEC beyond that of the strict Thomas- Fermi 
limit considered here. 

A. Local density perturbations 

Phonons have already been predicted 0, to cause 
instabilities in a homogeneous dipolar BEC when £dd > 
1. This can be seen directly from the Bogoliubov disper- 
sion relation between the energy Eb and momentum p 
for phonons in the gas 



EB = i{^) +^{l + £dd (3cos^-l)}^ 

(25) 

which can become imaginary when e^d > T indicating 
an instability. This dispersion relation has an angular 
dependence (9 is the angle between the momentum of the 
phonon and the external polarizing field) which further 
illustrates the richness of dipolar systems in comparison 
to the usual non-dipolar case. Equation (|25[1 is derived by 
adding to the Fourier transform, g, of the contact inter- 
action gS(r) that appears in the usual Bogoliubov disper- 
sion relation, the Fourier transform, Cdd{kikj — <5y/3), of 
the dipole-dipole interaction fllfl. This is an approxima- 
tion that assumes, as we do throughout this paper, that 
there is no screening of two dipoles by the other dipoles 
lying between them and also that the scattering of these 
two particles by the dipole-dipole interaction takes place 
within the Born approximation, as mentioned above in 

sec. nn 



In a trapped BEC with negligible kinetic energy, as 
considered here, we should expect an analogous insta- 
bility due to local density perturbations having a wave- 
length much smaller than the dimensions of the conden- 
sate. For example, Santos et al 0] recently showed that 
an infinite-pancake dipolar BEC, homogenous in two di- 
rections and parabolic in the third, is unstable when 
Edd > 1 for a density exceeding a critical value. 



B. Scaling deformations 

We use the term 'scaling' deformation to describe per- 
turbations that merely re-scale the parabolic solution 
1|11JI . i.e. that change R x = R y and R z from their equilib- 
rium values (|2(J[1 . but leave the basic form of the parabolic 
solution the same. Since the equilibrium values of the 
radii are determined by the transcendental equation 12111 , 
which also occurs in the context of a Gaussian variational 
solution, much of what we shall say below has already 
been described by other authors, including Santos et al 
|3, and Yi and You @- 

Information on the stability of the parabolic solution 
can be gained from analyzing the behavior of the energy 
functional E tot = E trap + E s - wavc + Edd evaluated over 
a general parabolic density profile (|llf> 

£"*m( 2 + 5) + t 1 - 

(26) 

in the vicinity of the solution (|20|l and 121|) and across 
the whole parameter space (k, Edd, 7)- One finds that for 
< Edd < 1 the solution given by the transcendental 
equation l|21|) is always stable, in the sense that it cor- 
responds to a global minimum of the energy functional. 
As Edd is increased and passes through unity, the solution 
matches smoothly onto one that is only metastable, i.e. it 
is only a local minimum in the energy landscape, and the 
global minimum is then a (prolate) collapsed state with 
k — > 0. Simultaneously with the turning of the stable 
into a metastable solution, a second branch of solutions 
appear at a smaller value of k, which correspond to un- 
stable saddle points that separate the local minimum of 
the metastable solution from the global minimum of the 
collapsed state at k = in the energy landscape. If one 
continues to increase Edd then one of two things happens, 
depending on the value of 7: if 7 is less than a critical 
value, 7 < 7 cr it = 5.1701, then, as £dd increases, even- 
tually the metastable and unstable solutions coalesce at 
£dd = e< dd> above which there are no solutions, not even 
metastable ones, and the energy landscape is just a con- 
tinuous slope down towards a collapsed state with k — 0. 
This critical value e^d* as a function of 7 is plotted in 
Figure If, however, 7 > 7 cr ;t then something rather 
surprising happens. As first remarked by Santos et al 
when 7 > 7 cr ;t there exists a solution metastable to 
scaling perturbations at a finite value of k for all values of 
Edd-, strictly speaking even for Edd = 00. (Note, however, 



CO 




FIG. 4: The critical value of the dipolar coupling, e^j , above 
which the condensate becomes strictly unstable — even the 
metastable state (local minimum in the energy landscape) 
no longer exists. However, above 7 cr i t = 5.1701, there is al- 
ways a solution metastable to scaling deformations even for 
arbitrarily large £dd (although not necessarily stable to local 
density perturbations — see text). 




p [units of R ] 

FIG. 5: Contour plot of the potential V (in units of the chemi- 
cal potential /_t) outside a condensate with « = 2 at £dd = 1.5. 
The closed contours drawn in dotted lines have V//i < 1; the 
lowest in the middle is at V / ji = 0.992 and the difference to 
the next higher is roughly 0.0017. The cut off contours out- 
side the condensate go from V/fi = 1.002 to 1.04 in steps of 
approximately 0.0034. 



that our value for 7 cr i t disagrees with that of Ref. p| but 
is that same as the one given in Ref. |j|). 

For completeness we would like to mention that prima 
facie the transcendental equation (|21|l for Edd > 1 has 
solutions also for n > 7. These come in pairs, one corre- 
sponding to a maximum and the other to a saddle point 
in the energy landscape. However, inspection reveals that 
these solutions have no physical relevance as for them the 
radius R x of the condensate, Eq. lj2U|) . comes out imagi- 
nary. 



C. Saturn-ring instability 

An examination of the dipolar potential outside the 
condensate, that is seen, for example by a single test 
atom placed beyond the boundary p 2 / 'R 2 + z 2 /R 2 = 1, 
reveals a new type of possible instability. Like the lo- 
cal density perturbations, it also does not preserve the 
parabolic form of the density profile. It turns out that 
for £dd > 1 the potential seen by atoms just outside the 
condensate exhibits a local minimum, i.e. the sum of 
trap and dipole-dipole potentials is locally lower than 
the chemical potential, which causes atoms to spill out 
from the condensate and fill this dip in the potential. 
Such an effect is peculiar to condensates with induced 
dipole-dipole interactions because these are long-range 
and thus give rise to a potential even outside the con- 
densate, whereas the potential due to s-wave scattering 
is short-range and thus zero outside the condensate. To 
investigate the dip in the outside potential we need to use 
the dipolar potential i|22l23|) that we calculated in Section 
Ivl As the expressions are rather awkward in Cartesian 
coordinates, we shall analyze them in spheroidal coordi- 



nates. 

The outside potential V is the sum of the trap potential 
and the dipole-dipole interaction potential $^° utl!ldc ) ( r ) _ 
The trap potential is positive and monotonically rising 
from the center. The dipole-dipole interaction potential 
outside the condensate is a solution of the (homogeneous) 
Laplace equation, i.e. <j> d ° utsldc ) ( r ) j g a harmonic function. 

Thus the Maximum Principle applies, and $|™ tsldc ) (r) 
must assume its maximum and minimum on the bound- 
aries of the domain, either at infinity or on the surface 
of the condensate. At infinity the dipole-dipole potential 
vanishes, which means that at large distances the total 
outside potential is positive and dominated by the trap 
potential. To ascertain whether the outer potential V has 
a local minimum one only needs to check whether the 
sum of trap and dipole-dipole potentials has a negative 
first derivative at the surface of the condensate in some 
outward direction. It is easy to see that local minima of 
V can occur only for 77 = 0: at a local minimum of V its 
first derivative with respect to £ and rj must vanish, but 
both the trap potential and the dipole-dipole potential 
are quadratic in 77, so that the first derivative dV/drj is 
proportional to 77 and thus vanishes only for rj = unless 
it vanishes for all rj. Therefore we only need to examine 
the derivative dV/d^ at 77 = 0; we find 



d_ 



2(1 -e d d) 



6?k 2 [1 - £dd/(«)] 



(27) 



where £b is the value of the spheroidal variable £ on the 
surface of the condensate, i.e. £b = l/Vl — k 2 for a pro- 
late BEC and £ B = 1/V« 2 - 1 for an oblate BEC. For 
oblate BECs the denominator in Eq. I|27(l is always pos- 
itive, and thus dV /d£ is negative if and only if Edd > 1- 
For prolate BECs the denominator is always positive in 



the stable and metastable regions of Fig |2 and thus in 
these regions dV/d^ is negative if and only if e dd > 1. 
Whenever <9V/<9£ is negative on the condensate surface, 
a local minimum of the potential lies somewhere outside 
the surface. The fact that this happens at r\ = means 
that the local dip in the potential occurs along a ring 
at z = around the condensate, and the flowing out of 
atoms from the main condensate into the dip causes the 
condensate to take on a Saturn-like appearance, which 
then leads to further instability. Figure gives an illus- 
tration of a typical potential, plotted as contours in the 
p — z plane. The flat part to the left of the center is the 
constant chemical potential inside the condensate. 




VII. ELECTROSTRICTION 



The volume V of the spheroidal BEC can be expressed 



as 



V = ^RlR z = ^ 
3 x 3 k 



Substituting R x from Eq. l]2()ll we can write it as 



V = 



4tt 



IbgN 

Annuo 2 , 



3/5 



,'3 « 2 /(«) 



(28) 



3/5 
(29) 



and then we can use the transcendental equation l|21(l to 
eliminate k in favor of 7 and £dd- In Fig. we plot V 
in units of [15gN/(4:TTmuJ x )] 3 ^ 5 as a function of Edd for 
various trap aspect ratios 7. The figure shows that con- 
densates in prolate traps and slightly oblate traps (with 
7 < 1.6630) get compressed by increasing dipolar interac- 
tions, while condensates in traps with aspect ratios above 
7crit = 5.1701 are being pulled apart as £dd rises. Be- 
tween 7 = 1.6630 and 7 cr i t is a range of trap aspect 
ratios for which the condensate is pulled apart at first 
but eventually compressed into collapse by higher values 
of the dipolar interaction strength. If, during an exper- 
iment, the condensate was imaged in the trap then the 
volume could be measured either directly from the radii, 
or by the central density which is inversely proportional 
to the volume, V = (5/2)N/n . 



VIII. RELEASE ENERGY 

If the trap is turned off the condensate expands ballis- 
tically and the s-wave and dipole-dipole interaction en- 
ergies are converted into kinetic energy, the so-called re- 
lease energy, which can be measured in an experiment. 
For the exact parabolic solution the release energy is 
given by 

E, cl = 15gN 2 [l - e dd f(K)]/(28nR 2 x R z ). (30) 

The release energy can also be expressed in terms of the 
chemical potential (|19f) as E xe ] = (2/7)7V/i, which is the 
standard expression for the Thomas- Fermi limit . 



FIG. 6: The volume V of the condensate for traps with aspect 
ratios 7={0. 1,0. 2, 0.5, 1,2, 3, 4, 5, 8} as a function of the dipole- 
dipole coupling strength Edd- At Edd = 0, i.e. without dipole- 
dipole interactions, the volume is proportional to 7 _2//5 , so 
that on the left edge of the graph higher volume corresponds 
to lower 7. That reverses with increasing e d d- 



IX. CONCLUSIONS 

The long-range anisotropic nature of dipole-dipole in- 
teractions gives condensates composed of dipolar atoms 
or molecules novel properties compared to those with 
only short-range interactions. Furthermore, the dipole- 
dipole interactions can be controlled in sign, magnitude 
and direction. We have shown that a simple inverted 
parabola remains an exact solution for the density-profile 
of a harmonically trapped dipolar BEC in the Thomas- 
Fermi limit by calculating the dipolar potential both in- 
side and outside the condensate region. The parabolic 
solution is stable, in the strict Thomas-Fermi limit, only 
for £dd < 1 ■ For £dd > 1 it is unstable against scaling per- 
turbations as seen in the energy functional, against per- 
turbations of different symmetry, such as phonons, and 
also due to a local minimum appearing in the potential 
outside the condensate. The effect of the dipole-dipole 
interactions upon the BEC is to change both its aspect 
ratio and its volume. The manner of this electrostriction 
depends on the aspect ratio of the external trap; as a 
function of the strength of the dipole-dipole interactions 
(relative to the s-wave scattering) the volume of the BEC 
can either increase or decrease, or even initially increase 
and then decrease at higher values of the coupling. 
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APPENDIX A: INTEGRATING OVER 
HOMOEOID SHELLS 

Here we demonstrate how the potential (|14f> can be 
obtained by integrating over successive thin shells, a 
method that works even in the general case of an el- 
lipsoid. The problem of the electrostatic/gravitational 
field within and outside a charged/massive ellipsoid is 
associated with some of the great names of nineteenth 
century mathematical physics. In his 1892 treatise on 
statics, Routh |26] notes that Chasles, Dirichlet, Jacobi, 
and Poisson all made contributions. Rodrigues is cred- 
ited with being the first to evaluate the potential of a 
solid ellipsoid of constant density, seemingly in 1815, and 
in 1833 Green [2]j "treated the subject in a very general 
manner" giving solutions for various cases of inhomoge- 
neous density, as in the current situation. Our derivation 
is adapted from Routh. Although we take a cylindrically 
symmetric density as input, the derivation holds for the 
general ellipsoidal case. 

Consider an ellipsoidal surface (x/R x ) 2 + (y/R y ) 2 + 
(z/R z ) 2 = 1, having semi-axes (R x , R y , R z ). A contin- 
uous family of concentric similar ellipsoids lying inside 
this outer surface can then be defined via the semi-axes 
(sR x ,sR y ,sR z ), where < s < 1. Note that by similar 
we mean that this family all share the same aspect ratios 
among their semi-axes, but are consequently not confo- 
cal. A surface which is confocal to the original surface 
obeys x 2 /(R 2 x + A) + y 2 /(R 2 y + A) + z 2 /(R 2 z + A) = 1, so 
that A = is obviously the original surface, and surfaces 
with A > lie outside the original. Considered in terms 
of the parameterization specified by s, the parabolic den- 
sity profile can be written as n = uq(1 — s 2 ). Since the 
cqui-dcnsity surfaces within the density profile are simi- 
lar, when computing the integral (|1U|I it makes sense to 
take as our basic volume element a thin "homoeoid" [2^ , 
defined as the shell bounded by two similar ellipsoidal 
surfaces, parameterized by s and s+ds, respectively. The 
volume of the thin homoeoid is dV = AirR x R y R z s 2 ds. 
When computing the potential at a point V(x,y,z) 
within a charged ellipsoid the contribution from the ho- 
moeoids interior to that point is of a different nature to 
that from those exterior to it, so we consider the two 
parts separately. 



1. Potential <f> m due to 'charge' interior to V 

Let dcf) denote the contribution to the total potential 
<f> from a thin homoeoid shell. In this section we require 
the potential d(f> m outside a thin homoeoid, labelled by 
its inner surface s, whose 'charge' density is n — no(l — 
s 2 ). The equi-potential surfaces outside a charged thin 
homoeoid s are confocal ellipsoids [26| 



V 



s 2 Rl + A 



s 2 R 2 y + A 



s 2 R 2 + A 



(Al) 



In our analogy (in which the dielectric constant eo = 
1) the 'potential' on a confocal surface A outside a thin 
homoeoid labelled by s, of density n and volume dV is 



n dV 
8tt 



du 



'A ^{s 2 R 2 x + u){s 2 R 2 y + u)(s 2 R 2 + u) ' 

(A2) 

This potential is most easily obtained by noticing that 
the 'charge' distribution on the homoeoid is identical to 
that of a solid ellipsoidal conductor with the same exter- 
nal dimensions and with the same total 'charge' (ndV) 
distributed over its surface. A homoeoid shell does not 
have uniform thickness, dh, being slightly thicker at the 
points furthest from the center 



dh 



s ds 




(A3) 



If, rather than a shell of thickness dh we consider the 
'charge' it contains to in fact be surface charge, of vari- 
able surface density uj(x,y,z), i.e. so that u) — ndh, we 
obtain exactly the same u(x,y,z) as in the well-known 
problem of the charged ellipsoid conductor (see, e.g. [28}), 
and thus the resulting potentials are also the same. 

Having established the potential due to a homoeoid 
shell we must integrate over all the shells s lying inside 
V, which is located on the confocal surface A. While V is 
of course fixed in space, the surface A passing through it 
is a different surface (i.e. different aspect ratio) for each 
homoeoid, i.e. A = A(s), becoming at the last instant 
a similar surface to the final homoeoid (upon which V 
itself sits). To simplify the algebra we put A = s 2 a in 
the equation for the ellipsoidal surface (|A1|) . We see that 
if s = then a = oo. We choose s = s to describe the 
similar ellipsoidal surface upon which V lies, and thus we 
have a = for s = s. To obtain the potential due to the 
charge interior to V we must evaluate <j) ln = J^ =Q d0' ; J 1 (s). 
Setting u = s 2 v in the integral l|A2(l . we find 



n R x RyRz 



f ds (l-s 2 )s f 

JO Jo 



(A4) 



dv 



M ^(R 2 x +v)(R 2 y +v)(R 2 +v) 
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which can be integrated by parts to give 
An _ n R x R y R z 



~2 ~4 
S 8 



dv 



o ^(R 2 +v)(R 2 J + v)(R 2 +v) 



2\Rl + a Rl + a R 2 + a 



A\Rl+(j Rl + <7 Rl + a 



da 

(R 2 x + a)(R 2 +a)(R 2 +a) 



(A5) 



where in the second term we have used Leibniz's formula 
for differentiating an integral and used Eq. IjAlfl to re- 
place s 2 and s 4 . 



Combining Eas . I A5I and I A6I we see that the parameter 
s drops from the sum 4> tot = cf> m + </> cx . In the general 
ellipsoidal case the remaining integrals in <f> tot can be ex- 
pressed in terms of elliptic integrals. In the spheroidal 
case (R x = R y ) they can be written in terms of more el- 
ementary functions. There are two cases to distinguish, 
depending upon whether the density profile is prolate (a 
cigar), R z > R x , or oblate (a pancake), R x > R z . If we 
define 



J{a,c) 



dv 



la (R 2 x +v)y/{Rl+v) 

7r 2 



V / v \ z • / 

i_ji+yi=Mm prolat£ 

-Rl \\-y/\-RHRl) 

arctan ( \J R 2 /R 2 — 1) oblate 



then the total potential at V is most compactly expressed 
as 



2. Potential OX due to 'charge' exterior to V 

The potential inside a homoeoid is much simpler to 
calculate since, just as in the case of the spherical shell, 
it is a constant, independent of position |29|. Returning 
to the solid ellipsoidal conductor model, the potential 
throughout a conductor is the same as on the surface 
defined by A = 0. Thus 



# ex (s) 



i dV 
~8^T 



du 



o .^(s 2 i?2 + u ){s 2 Rl + u){s 2 R 2 z + u) ' 



Note that because the lower limit of the integral is this 
time independent of s the ensuing treatment is simple. 
Integrating from the surface s, which includes V, out the 
boundary s = 1, we require 4> cx — Jg =s d<j) cx (s). Making 
once again the substitution u = s 2 v one immediately 
finds 



noR x R y R z 



i(l-.5 2 )-i(l-5 4 ) 
dv 



o ^(R 2 x + v)(R 2 J+ v)(R 2 z +v) 



(A6) 



uqR 2 R z 
2 

p 4 d 2 J 
8 d{Rlf 



dj 



4 + 2 d{R%) 
4 d 2 J 



dj 



z 



d(R 2 ) 

. , d 2 J \ 

3d(Rl) 2+pZ d(Rl)d{Rl)) 



(A7) 



and is identical to the Green's function result 114(1 . 
The relationship between the function 2 defined in the 
Green's function approach and J defined here is simply 
J = Z/R Z - 



APPENDIX B: THE POTENTIAL OUTSIDE THE 
CONDENSATE 



The calculation of the potential at a point V outside 
the boundary of the condensate follows very similar lines 
to that presented above for a point inside. The main 
difference is that one no longer has to consider 'charge' 
located exterior to V . The result is identical to Eq. (|A7|I 
except that the integral J is now given by 



J(a, c,A) = 



dv 



1 



A (R*+V)y/(R*+V) 



= < 



2 



In 



yxn?! + y/m, - ri 



y/^+Rl-y/Rl-Rl. 

arctan 



IRI-Rl 



X + R 2 



prolate 



oblate 



where A parameterizes an ellipse, confocal to the outer-boundary of the condensate, upon which the point V sits, and 
so obeys the quadratic equation p 2 /(R 2 + A) + z 2 /{R 2 + A) = 1. When solving for A one should take the positive 
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solution 

. _ p 2 - R\ + z 2 - Rl y 7 ( P 2 -Rl + z 2 - R 2 f + 4 ( P 2 R 2 Z + z2 fl 2 _ Rim 



Note that expressing the solution for the potential <f> in 
the form l|A7() requires that when taking the derivatives 
with respect to R x and R z then A is to be treated as a 
constant and is not to be differentiated. On the other 
hand, when going on to calculate the dipolar potential 
outside the BEC, $dd = — Cdd e»ejViVj<£(r), then A is 
to be treated as a function of (p, z), and should be differ- 
entiated. This makes the exact expression for $dd out- 
side the condensate complicated. However, in the limit 



p 3> R x , z 3> R z , one may use the asymptotic result 



outside N (R 2 x -R 2 z )(p 2 -2z 2 ) + U(p 



8 77 



7{p 2 + Z 2 ) 5 /2 



(Bl) 



which turns out in practice to be remarkably accurate, 
even right up to the condensate surface providing the 
condensate is not too aspherical (in which case the next 
terms in the multipole expansion should included.) 
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